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The problem of Josephson current through Coulomb-blocked nanoscale superconductor-normal- 
superconductor structure with tunnel contacts is reconsidered. Two different contributions to the phase-biased 
supercurrent I(ip) are identified, which are dominant in the limits of weak and strong Coulomb interaction. 
Full expression for the free energy valid at arbitrary Coulomb strength is found. The current derived from 
this free energy interpolates between known results for weak and strong Coulomb interaction as phase bias 
changes from to n. In the broad range of Coulomb strength the current-phase relation is substantially 
non-sinusoidal and qualitatively different from the case of semi-ballistic SNS junctions. Coulomb interaction 
leads to appearance of a local minimum in the current at some intermediate value of phase difference applied 
to the junction. 

PACS: 73.21.-b, 74.50.+r, 74.45.+C 
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The Josephson current in the contact of two bulk 
superconductors through a small normal grain (SINIS 
structure) was recently studied in £Q within the saddle- 
point approximation for the effective action functional 
which described [5] superconductive proximity effect in 
the presence of Coulomb interaction. We found in £Q 
that Coulomb blockade in the grain results in a very 
special Josephson current dependence on the phase dif- 
ference ip (see dashed lines in Fig.0. As ip approaches it 
the proximity effect in the grain is suppressed, Coulomb 
blockade becomes relatively strong, and spectral mini- 
gap induced in the normal grain becomes exponentially 
small. This lead us to erroneous conclusion that su- 
percurrent through the grain is exponentially weak as 
well. In fact, as was pointed out to us by I. S. Be- 
loborodov and A.V. Lopatin [3], supercurrent in the 
SINIS structure may flow without any spectral mini- 
gap at all. The same result was obtained originally by 
C. Bruder, R. Fazio and G. Schon in 0] by means of 
perturbative analysis that is valid as long as charging 
energy Ec = e 2 /2C is much larger than the proximity- 
induced minigap (denoted below as E g ). This additional 
(with respect to our result in JT|) contribution to the su- 
percurrent is due to nuctuational diffuson and Cooperon 
modes in the N grain, as explained below. 

In pp the replicated dynamical sigma-model was used. 
The Coulomb interaction in the grain is taken into ac- 
count in the framework of the adiabatic approximation 
developed in . The key point of this approximation is 
the separation of energy scales: the electric potential of 
the grain fluctuates at frequency much larger than the 
proximity induced minigap. This assumption is valid 
provided Ec 3> 5, where 5 is average level spacing in 
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the grain per one spin projection. The saddle point 
of the sigma-model gives the free energy of the system 
Fq(lp); the current is than calculated using the identity 
I = (2e/h)dF /dip. 

The above-mentioned additional contribution to the 
supercurrent is due to the fluctuations near the saddle 
point of the sigma-model. Below we present a somewhat 
more general result calculating nuctuational correction 
to the total free energy of the system. The supercurrent 
being the derivative of the free energy with respect to 
ip acquires significant correction in the regime when the 
saddle point itself gives exponentially small result. This, 
anyhow, happens when the phase difference ip comes 
close to 7t. Therefore, the results of pQ and for the 
Josephson current are, in fact, valid in two limiting cases 
of weak and strong Coulomb effect correspondingly. 

In order to find Josephson current in the full range of 
Coulomb/proximity ratio, it is necessary to supplement 
our results presented in £Q by nuctuational contribution. 
The saddle-point approximation used in m2 is justified 
by the inequality E g 3> 8. The correction due to fluc- 
tuations near the saddle point is negligible in this limit. 
When the phase bias ip is close to it the parameter E g 
becomes exponentially small and the above inequality 
is violated. It is this violation that makes nuctuational 
correction important. However, due to the large value of 
junction's dimensionless conductance, it is sufficient to 
consider the nuctuational contribution in the Gaussian 
approximation not going beyond the quadratic expan- 
sion of the action in soft modes. In this paper we extend 
the approach of ^ E] to allow for Gaussian fluctuations 
near the sigma-model saddle point. The fluctuations are 
calculated on the background of a non-zero proximity- 
induced minigap. The resulting dependence I(ip) inter- 
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polates between that from pQ and from 0] as cp varies 
from to 7T. Moreover, in the crossover region a local 
minimum of supercurrent appears (see Fig. 

We start with the derivation of the fluctuational con- 
tribution to the free energy valid for arbitrary rela- 
tive strength of Coulomb blockade and proximity effect. 
Then the current-phase dependence is found by numeri- 
cal differentiating. To reduce unnecessary complications 
the temperature is put to zero. We assume the SINIS 
junction between two bulk superconductors with tunnel 
contacts characterized by large (in units e 2 /K) normal 
conductances Gl and Gr. It is convenient to introduce 
the effective conductance 

G( V ) = tJg 2 l +G 2 r + 2G l G r cos^. (1) 

Formally, one may treat the system as an SIN junction 
with one superconductive lead and normal conductance 
given by the above expression . We quantify the prox- 
imity effect in the normal grain by the bare value of 
induced minigap E g (ip) — G((p)5/A. This minigap is 
realized if the Coulomb interaction is absent. 

The sigma-model for SINIS junction deals with 
the matrix field Q^, that bears two replica and two 
Matsubara energy (or, equivalently, imaginary time) 
indices along with particle-hole structure in Nambu- 
Gor'kov space. Another field is the scalar phase 
dependent on the imaginary time and replica number. 
The action has the form 



pl/T 

S[Q, K] = -- Tr( £ f 3 Q) dT 



4E C 



- G{ip) tr |Q~ (fi cos 2K% + f 2 sin 2K?) 



(2) 



The symbol fi is used for Pauli matrices operating in 
Nambu-Gor'kov space. The operator 'Tr' implies sum- 
mation over all indices including Matsubara energies 
and replicas, while 'tr' denotes trace in Nambu-Gor'kov 
space only. 

The adiabatic approximation allows to integrate out 
the field K assuming Q to be fixed. The steady matrix 
Q is diagonal in energies and trivial in replicas 



2ir6 ab 5(e - e') 



£T 3 + Eafl 



+ E 2 



(3) 



Here E g is the proximity-induced minigap in the normal 
grain renormalized by the Coulomb interaction. The 
value of Eg will be determined self-consistently later. 



Using the adiabatic approximation we derive the ef- 
fective Hamiltonian that determines the dynamics of 
phase K 



H = E C [-d 2 /dK 2 - 2qcos2K] , 
^_E g {^)E g 



Jg 2A 

i^r l0g i7 



(4) 
(5) 



The parameter q quantifies the relative strength of 
Coulomb blockade in comparison with proximity ef- 
fect. In the limit g > 1 Coulomb interaction is effec- 
tively weak and can be treated perturbatively, wile in 
the opposite case q -C 1 Coulomb blockade destroys the 
proximity effect up to an exponentially small correction. 
The value of q is determined self-consistently along with 
E g . 

At zero temperature K is frozen in the ground state 
of the Hamiltonian (jJJ with energy Eo(q). Then the 
total free energy acquires the form 



Fo(<p) 



e de 



EM- 



(6) 



The divergent integral is to be regularized by subtract- 
ing its value for the "normal" state with E g = 0. In all 
subsequent analysis we assume this regularization to be 
done. 

The minigap E g is set by the condition dFo/dE g = 0. 
This gives the self-consistency equation 



E a 



1 dE 

E g (cp) 2E C dq 



(0|cos2i<:|0). 



(7) 



Last expression implies average value at the ground 
state of (@J. Together with © this equation forms a 
closed system that determines q and E g . 

In pQ we calculated the supercurrent using the iden- 
tity Iq = (2e/h)dF /d(p. The result was 



E,. 



e5_f 

m\E 6 



2A 



G L G R sin <p log 



E„ 



(8) 



This dependence of supercurrent on <p is shown in Fig. ^ 
by dashed lines. 

To take into account fluctuations of Q near the found 
saddle point (for detailed calculation see [5]) we use 
the parametrization of Q as rotated fi matrix: Q = 

y-l e -iW/2f ie iW/2y_ Thig chok;e j g motivated by thc 

fact that Q — fi at e — 0. Below we'll see that 
different modes of fluctuations near the saddle point, 
diffusons and Cooperons, decouple in this representa- 
tion. The diagonal in energies and replicas matrix V 
is determined by the identity Q = V~ x f\V and ex- 
pressed as V = cos(7r/4 — 9/2) — if2sin(7r/4 — 9/2) 
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with 8 being a standard Usadel angle, tanf? = E g /e, 
dependent on Matsubara energy. The matrix W de- 
scribes deviations of Q from the saddle point Q. W 
anti-commutes with f\ and hence contains two compo- 
nents W^, = T 3 df £ , + T 2 cf £l . Below the (off-) diagonal 
in Nambu space element d ££ , (c ££ ,) is referred to as dif- 
fuson (Cooperon) mode. However, these modes are only 
analogs of standard diffuson and Cooperon that describe 
fluctuations near normal metallic saddle point with no 
minigap. 

Now we substitute the above parametrization into l|2|) 
and expand the action to the second order in W. The 
result is a sum of three terms 

S[W, K] = So [K] + S x [W, K] + S 2 [W, K] , (9) 

where So [K] is the action corresponding to the Hamilto- 
nian (0J. The terms Si,2[W, K] are linear and quadratic 
in W respectively; the explicit expressions for them can 
be found in 0. Our strategy is to integrate out the 
phase K. As the fluctuations near the saddle point are 
assumed to be small, we treat the last two terms of © 
as a perturbation to the bare action So [if] or, equiv- 
alently, to the Hamiltonian Q. Thus we expand the 
statistical weight e~ s ^- W ' K ^ to the second order in Si 
and to the first order in S 2 - Integration with respect to 
K implies averaging of these terms at the ground state 
of Q. Once the integral is calculated we rewrite the 
result in the form of a single exponent 



DWDKe- s ^=e~^ / DW e 



j(l)_Q< 2 >_C< 2 > 



(10) 

s« = (s 1 ), sf = (s 2 ), s% = - {Sx)2 . 

(11) 

The integral of the term e~ So yields a W-independent 
factor to the partition function corresponding to the free 
energy F calculated at the saddle point. The value of 
Fq is given by © while N is the number of replicas. 
The symbol (...) denotes the average with respect to 
the ground state of the Hamiltonian To calculate 
the terms S^ and Sq 2 ' we use the identity (sin 2K) = 
while the average value of cos2if is determined by (7J. 
The self-consistency equation provides S^ = as it 
should be at the saddle point. The two remaining terms 
describe fluctuations near the saddle point. They can 
be written in the form 



S 



(2) 
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dede' 
8tt<5 



E l 



ab ba 
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,a e'e 



(12) 



c(2) _ 
°int — 



G 2 y—^ f de de' dto 

x [\ c {e 1 e'-oj)cZ l c7 + ^ 

+ \ d (e,e';u)d a e %,d a e ? +u , e+u 



(13) 



In the last expression we use the following notations 



A c (e, e'; uj) — — X(e — e') cos • 



cos ■ 



(14) 



— Y(e — e ) cos 
4 ^ 



(15) 

X M = E l(0lcos2X|,)| ^ 2 2 ^-_^ (16) 

n>0 

yM^^|(0|sin2if|n)| 2 2 ^"' g ° } (17) 

The two functions X(ui) and Y(u>) appeared from the 
averaging of S 2 . They are nothing but the Fourier com- 
ponents of irreducible correlators ((cos 2Kq cos 2K T )) 
and ((sin2ifosin2.KV)) respectively. In the expres- 
sions l|T5fr and (|T7J) we use \i) and Ei to denote eigen- 
vectors and eigenvalues of JJJ. 

With all these definitions in hand we employ the stan- 
dard replica trick to calculate the free energy 



F = F +T lim — 

n^o N 



DWe~ b <> 



(18) 



The last term of this expression is the fluctuational con- 
tribution. It contains Gaussian integral with rather 
complicated quadratic form in the exponent. The term 

(2) 

S is fully diagonal with respect to both replica and 
Matsubara energy indices of c and d components. The 

(2) 

complication arises from the S int term where differ- 
ent energies are coupled (|13fl . The Gaussian integra- 
tion is equivalent to computing the determinant of this 
quadratic form. To find the value of this determinant we 
use the standard trick |Hj . Let us consider the derivative 
of free energy with respect to G 2 : dF/d(G 2 ). The fac- 
tor G 2 is present explicitly only in the term. Then 
the derivative is 



dF 



= lim 



T 



d(G 2 ) n^o NG 2 



DWS$e- s o 



(19) 



The pre-exponent of the integrand is quadratic in c and 
d. Hence the differentiating with respect to G 2 reduced 
the problem of calculating the determinant to the calcu- 
lation of the inverse matrix elements. They correspond 
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to the two-particle Green functions, that are Cooperon 
and diffuson. The Cooperon is 



The analogous screening function for the correlator 
D ab (e 7 e';uj) is 



(c ab ,(?? 



■uj ,e+uj 



)= .. 

= 2ir6 aq 



e'b- 



-e'+ui'p 
■ e+uj, q 
ib 



uj-uj')C ab {e,e'-uj). (20) 



The analogous diffuson function is similar to the above 
average with c components being replaced by d. We 
denote this propagator by D ab (e, e'; uj). Angle brack- 
ets in the last expression imply the average with the 

(2) (2) 

Gibbs weight given by the quadratic action Sq ' + S^ n ( . 
This averaging also contains normalization by the parti- 
tion function that is the determinant we are calculating. 
However, this normalization factor is canceled when the 
N — > limit is taken in (|19() . (Note, the replica trick 
was originally invented namely for this cancelation.) 

Below we concentrate on the quantity C ab (e, e'; uj). 
Another propagator is found in analogous way; the 
only change is the replacement of A c by A^- To find 
C ab (e, e'\ uj) we have to solve a simple Dyson equation 



(21) 



where the bare correlator Co (£,£') is determined by in- 

(2) 

verse eigenvalues of the diagonal quadratic form Sq 
and the vertex is the matrix element of S: 



(2) 



C (e,e') 



e'.b 



e 2 + El 



El 



(22) 



V =G 2 S ab X c (s,e';uj). (23) 



The solution to the equation l|21|) is 



C ab {e,e'-uj) = C {e,e') 



2ttS(uj) 



™u G^ 1 \ c (e. £■': uj) „ , . . 
6 l-Jc{e-e>) C0ie+UJ > E + U) . 



(24) 



In the denominator of the last term the screening func- 
tion C(2Q) appears 



C(20) 



de 
2^ 



A c (e + O, £ - Q; 0)C o (e +Q,e-Ct) 



-X{2fl) 



2A E 2 O 

log -= \ 1 + — f arcsmh — 

E n V ^ 2 E„ 



(25) 



This integral contains logarithmically divergent contri- 
bution that gives the first term in square brackets. The 
second term comes from small values £ < max{ri,_B g }. 



V{2Q) = -y(2fi) 



, 2A arcsinh(fi/S a ) 
log — 



E„ 



i + E 2 /n 2 



(26) 



Now everything is ready for the calculation of the 
derivative dF/d(G 2 ). The pre-exponent in (|19|l con- 
tains the sum over single replica index. The saddle point 
is trivial in replica space. Therefore, this sum will be 
canceled by 1/N. Another feature of pre-exponent is 
the factor 2irS(0). This factor appears from the delta- 
function in the definition of Cooperon (|20|) . The same 
is also true for the diffuson term in the pre-exponent. 
At finite temperature 2tt5(Q) = l/T that cancels tem- 
perature in (|19|l . Once this cancelation is established 
we can safely assume T = 0. 

After substitution of (|24|l and similar expression for 
D ab (e, e'; ui) into (|19|l and integration with respect to uj 
and to the sum e + e' we come to a single integral over 
Q = £ - e' 



dF 
9(G2) 



dn 

47T 



C(Q) 



V(Q) 



1 - G 2 C(fl) 1 - G 2 V(n) 



(27) 

In the limit G 2 = the vertex part of the action, sj^. , is 
absent. This leads to the absence of fluctuational con- 
tribution at G 2 = in the limit N — > 0. Using this fact 
we finally come to the full expression for the free energy 
by integrating l(2*7|l 



F = F o + J lo g[! - G 2 c(n)] [i - G 2 v(n)] . 



(28) 



The total Josephson current is now easy to find by 
differentiating free energy I = (2e/h)dF/dip. The cur- 
rent derived from the first term Fq was found in £Q. 
Rather simple expression (JSJ exists for this quantity. 
The fluctuational contribution is much more compli- 
cated. The dependence on phase difference ip is con- 
tained not only in the factor G 2 according to JU but 
also in the screening functions. The situation is very 
much simpler in the physically interesting limit of strong 
Coulomb blockade. The minigap is strongly suppressed 
and the saddle point itself produces negligible contri- 
bution to Josephson current. In the fluctuational part 
of the free energy we put E g — 0. The expression J3J 
becomes a free particle Hamiltonian as q — 0. Only 
first term is left in both sums (|16|l and Q17[). hence 
X(uj) = Y(uj) = AE C /(IQE 2 C + uj 2 ). The screening 
functions (|25|l and (|26|l become identical and contain 
only log(A/f2) in square brackets. The dependence of 
free energy on ip is now provided only by the factor G 2 
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Fig 1. The dependence of Josephson current / on phase 
bias ip. The current is normalized by its maximal value 
that is reached at ip = n/2 in the absence of Coulomb in- 
teraction J max = {eS/4h)G L G R log(8A/<VG| + G%). 
Dashed lines show the dependence JSJ found in 
without fluctuational correction. This correction is de- 
picted by the dotted lines, while solid curves are sums 
of the two contributions. The three plots correspond 
to different relative strength of Coulomb blockade: (a) 
weak interaction Ec6/Eg(0) — 0.5, (b) intermediate 
E c 5/E 2 g (0) = 1.5, (c) strong blockade E c 8/E 2 g (0) = 
2.5. For all three plots we assume symmetric junction 
with Gl = Gr = 20 and A/5 = 3000. As ip approaches 
7r the Coulomb interaction always becomes strong and 
the current is dominated by the fluctuational contribu- 
tion fM . 



in H28fl . For the Josephson current in strong Coulomb 
blockade regime we have 

Hf) = ^ G L G R sin <p log (29) 

Thus the result of 0] is reproduced. 

In the opposite limit of weak Coulomb blockade the 
fluctuational contribution to the supercurrent is small 
in comparison with Iq. When the phase difference 
changes from to 7r the system goes from weak to 
strong Coulomb blockade regime. This means that the 
result © gradually transforms into (|29[1 . Numerical 
differentiating of H28|l gives the solid curves plotted in 
Fig-fflfor the current-phase dependence. Note the local 
minimum that appears in the crossover region. There is 
no simple analytic theory for this effect. However, this 
is likely to be the most prominent feature of the system. 

In conclusion, we have calculated the fluctuational 
correction to the free energy and to the Josephson cur- 
rent in Coulomb blockaded SINIS junction. This cor- 
rection plays major role in the limit of strong Coulomb 
interaction. At intermediate value of phase bias a well- 
defined local minimum of Josephson current appears. 
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the contract RI-112/001/417, and by RFBR grant 04- 
02-16348-a. P.M.O acknowledges the support of "Dy- 
nasty Foundation" . 
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